Estatística Avançada - Projeto Final - FGV

CO2 Na França

Grupo:

• Daniel Mello Duarte Morais
• David da Guia Carvalho
• Elaine Maria Lopes Loureiro
• Guilherme Vieira Dantas
• Mallena Ferreira de Morais Costa

Índice

Resumo

No presente projeto, iremos, a partir de base de dados disponibilizada em aula, gerar um modelo linear para a previsão de CO2 em KT. Os resultados obtidos foram positivos e seguiram todos os requisitos necessários para a validade da regressão:

• Um p-valor acima de \(5\%\) para os testes de Anderson-Darling e Shapiro, o que indica normalidade em nossos resíduos, conforme o esperado.
• Um p-valor também elevado no teste de Breusch-Pagan, o que nos permite aceitar a hipótese de Homocedasticidade ainda que a reta do gráfico Scale-Location não seja propriamente horizontal.
• Um p-valor elevado para o teste de Durblin-Watson.
• Tudo isso a um \(R^2\) na ordem de \(99\%\)

O modelo foi desenvolvido, no âmbito de nosso grupo, para os dados relativos à França. Nas próximas seções, iremos descrever cada um dos passos realizados, desde a análise exploratória de dados até os resultados finais, em formato de storytelling.

1. Primeiro Passo - Verificação de Dados

Iremos iniciar o trabalho verificando os dados que temos para a realização da análise. Para isso, seguiremos a seguinte sequência de passos:

• Verificar o arquivo de entrada
• Conferir quais atributos possuem informações para todos os instantes de tempo da análise
• Conferir quais atributos possuem informações parciais (missing data)
• Conferir e eliminar atributos que não possuem informações ou possuem uma quantidade muito reduzida de dados, de forma a inviabilizar a análise

setwd('~/Área de Trabalho/stu/MBA - FGV/Estatística Avançada/Trabalhos/Trabalho-Final')
df_in <- read_excel('./FRA_Country_en_excel_v2.xls')

## New names:
## * `` -> ...3
## * `` -> ...4
## * `` -> ...5
## * `` -> ...6
## * `` -> ...7
## * … and 53 more problems

df_in[, 1:3] %>% head

Os dados parecem extremamente desorganizados. Precisamos, primeiramente, eliminar as duas primeiras linhas e tornar a terceira linha o header do Data Frame:

colnames(df_in) <- df_in[3,] %>% unlist
df_in <- df_in[-c(1, 2, 3),]
df_in[, 1:3] %>% head

Sabemos que o país é a França, então podemos eliminar o nome e o código do país. O nome do indicador será utilizado para se referir à variável e o nome do indicador será eliminado da tabela. Um segundo Data Frame relacionando os códigos dos indicadores aos nomes será criado e funcionará como uma tabela de metadados:

df_in_meta <- df_in %>% select(`Indicator Name`, `Indicator Code`)
df_in_meta %>% head

Eliminando os metadados e as informações redundantes de país da tabela de dados:

df_in <- df_in %>% select(-`Country Name`, -`Country Code`, -`Indicator Name`)
df_in[, 1:3] %>% head

Podemos notar que existem muitos indicadores com um altíssimo número de dados perdidos. Iremos então verificar em uma tabela a quantidade de dados perdidos por indicador, para descobrir quais indicadores não serão levados em consideração na análise:

df_in$`N. Missing` <- apply(df_in, 1, function(x) (is.na(x) %>% sum))
df_in$`Perc. Missing` <- df_in$`N. Missing` / (ncol(df_in) - 2)

df_in_missing <- df_in %>% 
  select(`Indicator Code`, `N. Missing`, `Perc. Missing`) %>% 
  arrange(desc(`N. Missing`))

Escrevendo a a saída em uma tabela Shiny formatada e em um gráfico de barra com o percentual de dados perdidos:

ggplot(data=df_in_missing, aes(x=`N. Missing`)) + 
geom_histogram(color='darkgreen', fill='white') + theme_minimal() +
  labs(x='Quantidade de dados perdidos',
       y='Freq. Absoluta',
       title='Missing Data',
       subtitle='Analisando quantidade de dados perdidos por tipo de variável')

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

df_in_missing

Se tomarmos um limiar igual a 5, admitiremos um valor de perda máximo de dados igual a a aproximadamente 10% e teremos no mínimo 50 pontos de amostra para gerar o modelo:

ggplot(data=df_in_missing %>% filter(`N. Missing` <= 5), aes(x=`N. Missing`)) + 
geom_histogram(color='darkgreen', fill='white') + theme_minimal() +
  labs(x='Quantidade de dados perdidos',
       y='Freq. Absoluta',
       title='Missing Data',
       subtitle='Filtrando gráfico anterior para regiões com menos de 6 dados perdidos')

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

df_in_missing %>% filter(`N. Missing` <= 5)

De acordo com o histograma, ainda que tal critério rigoroso seja adotado, uma grande quantidade de indicadores continua disponível: reduzimos o escopo da análise de 1200 indicadores para 200 indicadores e, ainda assim, temos um grande volume de dados.

Podemos, então, plotar conjuntamente os diagramas para verificar, de maneira superficial, o aspecto das curvas que serão submetidas a análise. Para facilitar tal tarefa, iremos transpor a tabela filtrada:

df_in <- df_in %>% filter(`N. Missing` <= 5)
data_matrix <- df_in[1:nrow(df_in), 2:(ncol(df_in) - 2)]
df_in_t <- t(data_matrix) %>% as.data.frame
colnames(df_in_t) <- df_in$`Indicator Code`
df_in_t$Year <- colnames(df_in)[2:(ncol(df_in) - 2)]
df_in_t[,1:3] %>% head

Iremos eliminar o ano de \(2015\) da análise pois foi observado que a coluna de tal ano não possui nenhum valor para nenhuma variável.

Plotando a evolução de cada uma das variáveis em função do tempo:

df_in_p <- df_in_t %>% gather('Indicator', 'Value', -Year)

## Warning: attributes are not identical across measure variables;
## they will be dropped

df_in_p <- df_in_p %>% filter(Year != '2015')

df_in_p <- df_in_p %>% transform(Year=as.numeric(Year)) %>% transform(Value=as.numeric(Value))
ggplot(df_in_p, aes(x=Year, y=Value, color=Indicator)) + geom_line() + 
  theme(legend.position='none', panel.background = element_blank(),
        panel.grid.major = element_line(colour = 'gray'), 
        panel.grid.minor = element_line(colour = 'gray')) +
  labs(x='Ano',
       y='Valor',
       title='Evolução das Variáveis Escolhidas',
       subtitle='Uma primeira análise das séries temporais')

## Warning: Removed 148 rows containing missing values (geom_path).

Podemos observar que existem dados perdidos devidos aos “Warnings” recebidos. Na próxima seção iremos corrigir esse problema.

2. Segundo Passo - Preenchimento de Missing Data

A função “approx” da biblioteca zoo permite que realizemos interpolação linear para preencher os dados faltantes de cada uma das séries. Esse procedimento é mais coerente que simplesmente copiar os dados do passado ou do futuro para completar as células iguais a NA.

Esse procedimento não considera os pontos extremos com valores iguais a NA, razão pela qual, após a interpolação, escrevemos uma rotina para preencher os valores extremos com o próximo valor nao nulo ou o valor não nulo imediatamente anterior.

indicators_list <- colnames(df_in_t %>% select(-Year))
for (curr_indicator in indicators_list) {
  
  curr_list <- df_in_t[curr_indicator]
  curr_list <- na.approx(curr_list, na.rm = FALSE)
  
  non_NA_index <- which(!is.na(curr_list))
  first_non_NA <- min(non_NA_index)
  last_non_NA <- max(non_NA_index)
  
  if (first_non_NA > 1){
    curr_list[1:(first_non_NA - 1)] = curr_list[first_non_NA]
  }
  if (last_non_NA < length(curr_list)) {
    curr_list[(last_non_NA + 1):length(curr_list)] = curr_list[last_non_NA]
  }
  
  df_in_t[curr_indicator] = curr_list
}

Finalmente, podemos plotar novamente os gráficos:

df_in_p <- df_in_t %>% gather('Indicator', 'Value', -Year)

## Warning: attributes are not identical across measure variables;
## they will be dropped

df_in_p[is.na(df_in_p)] <- 0

df_in_p <- df_in_p %>% transform(Year=as.numeric(Year)) %>% transform(Value=as.numeric(Value))
ggplotly(ggplot(df_in_p, aes(x=Year, y=Value, color=Indicator)) + geom_line() + 
  theme(legend.position='none', panel.background = element_blank(),
        panel.grid.major = element_line(colour = 'gray'), 
        panel.grid.minor = element_line(colour = 'gray')) +
  labs(x='Ano (Variáveis reativas: encostar mouse para verificar)',
       y='Valor',
       title='Variáveis Corrigidas'))

Não há alterações visíveis pois os dados perdidos se concentraram majoritariamente em posições extremas das seŕies. Em todo caso, os warnings de dados faltantes não se encontram mais presentes e podemos prosseguir para a próxima etapa: analisar a correlação e o relacionamento entre as variáveis.

3. Terceiro Passo - Análise de Correlação

• Em um primeiro momento, estudaremos os relacionamentos mútuos entre cada par de variáveis que podemos tomar.
• Em seguida, iremos verificar as correlações de cada uma das variáveis com a variável que iremos estudar (a emissão de CO²).
• Finalmente, com essas informaçoes, iremos organizar e categorizar grupos de variáveis que estão fortemente correlacionadas utilizando técnicas de clusterização.

Etapas da Análise de Correlação:

3.1. Correlações Mútuas entre Variáveis

Podemos plotar a correlação entre as variáveis em um mapa de calor. É inviável escrever o nome de cada uma das variáveis na matriz, entretanto, podemos, ao menos, observar o aspecto de tal mapa.

cor_mat <- df_in_t %>% select(-Year) %>% cor
ggcorrplot(cor_mat, tl.cex=0) +
  labs(title='Matriz de correlação como Mapa de Calor',
       subtitle='Identificando se há multicolinearidade')

Aparentemente, os grupos de variáveis são fortemente relacionados entre si. Isso indica que podemos agrupar as variáveis em um número muito mais compacto de categorias sem que isso signifique prejuízo em nossa análise.

Vejamos a correlação entre cada uma das variáveis com o dado em estudo na próxima subseção.

3.2. Correlações com Emissão de CO²

Antes de agruparmos as variáveis em categorias correlacionadas, vamos dar um “Zoom” na análise da seção anterior e nos ater à variável em estudo. Trata-se da emissão de CO² em KT, representada pela variável EN.ATM.CO2E.KT.

Podemos visualizar as correlações de cada uma das demais variáveis por meio de um gráfico de barra.

x_plot <- rownames(cor_mat)
y_plot <- cor_mat[, 'EN.ATM.CO2E.KT']
df_plot <- data.frame(X=x_plot, Y=y_plot)

ggplotly(ggplot(data=df_plot, aes(x=reorder(x_plot, -abs(y_plot)), y=abs(y_plot))) +
  geom_bar(stat='identity', color='darkgreen', fill='white') + 
  theme_minimal() + theme(axis.text.x=element_blank()) +
  labs(x='Nome da Variável (Gráfico reativo: encostar o mouse para verificar)', 
       y='Correlação com Emissão de CO2 (KT)',
       title='Correlações com Variável Alvo'))

Plotando o valor absoluto da correlação:

ggplotly(ggplot(data=df_plot, aes(x=reorder(x_plot, -abs(y_plot)), y=abs(y_plot))) +
  geom_bar(stat='identity', color='darkgreen', fill='white') + 
  theme_minimal() + theme(axis.text.x=element_blank()) +
  labs(x='Nome da Variável (Gráfico reativo: encostar o mouse para verificar)', 
       y='Correlação com Emissão de CO2 (KT)',
       title='Correlações com Variável Alvo'))

Não é, a princípio, necessário saber o nome de cada uma das variáveis representadas por cada uma das barras. O que sabemos, por hora, é que a correlação varia praticamente de forma contínua e todo tipo de correlação ocorre no grande conjunto de dados que temos disponíveis.

Assim, na próxima seção, iremos tentar agrupar as variáveis em grupos correlacionados no intuito de:

• Explicar quais variáveis possuem relacionamentos mútuos
• Utilizar tais relacionamentos para criar um modelo explicado por um menor número de variáveis explicativas.

Antes de efetuarmos o agrupamento, vamos retirar de nossa tabela a variável que queremos explicar>

df_out <- df_in_t[, 'EN.ATM.CO2E.KT']
df_in_t <- df_in_t %>% select(-EN.ATM.CO2E.KT)

3.3. Agrupamento de Variáveis Correlacionadas

As técnicas de aprendizado não supervisionado consistem na categorização de variáveis por meio de algoritmos de Cluster. Esse tipo de procedimento pode ser aplicado no agrupamento de variáveis fortemente correlacionadas.

Para isso, podemos imaginar que quanto mais correlacionadas duas variáveis forem, mais próximas elas estarão no hiperespaço de variáveis explicativas. Ou seja: podemos adotar uma métrica \(\mathcal{F}(Cor(X, Y)) = \mathcal{D}(X, Y)\) que corresponde à distância entre os vetores \(X\) e \(Y\).

Essa métrica é tal que, quanto mais correlacionadas forem as variáveis, menor será a distância entre elas, de tal forma que variáveis \(100\%\) correlacionadas serão separadas por uma distância igual a zero.

Além disso, a correlação será avaliada em valor absoluto neste ponto. Isso ocorre porque variáveis com correlações próximas a \(-1\) também podem ser consideradas extremamente próximas, ainda que as variações de cada uma delas ocorram com sinais trocados.

Definiremos a função em questão como: \(\mathcal{F}(Cor(X, Y)) = \mathcal{D}(X, Y) = 1 - |Cor(X, Y)|\).

A correlação entre cada par de variáveis formará uma matriz de ordem \(N \times N\), onde \(N\) é a quantidade de variáveis presentes na base de dados (em torno de 200 variáveis).

Assim, essa matriz terá a forma: \(\mathcal{D}_K(X_i, X_j) = 1 - |Cor(X_i, X_j)| = (d)_{ij} = D_{N \times N}\)

Os pontos serão agrupados por meio da técnica de clusterização hierárquica, considerando-se a matriz de distância definida acima: \(\mathcal{F_K}(Cor(X, Y)) = \mathcal{D}(X, Y) = 1 - |Cor(X, Y)|\).

A clusterização hierárquica não pede que o usuário determine o número de clusters desejado. Ao contrário, o usuário pode visualizar um dendograma com diferentes níveis de agrupamento e escolher aquele que melhor se aplica à necessidade do problema:

hclust_obj <- (1 - abs(cor_mat)) %>% as.dist %>% hclust
dend <- hclust_obj %>% as.dendrogram
dend %>% dendextend::set('labels_color', 'white') %>% plot(xaxt='n')

Precisamos escolher uma altura na qual a árvore será cortada. Para facilitar essa análise, plotemos o número de clusters em função da altura de corte:

height <- seq(from=0, to=1, by=0.001)
n_clusters <- sapply(height, function(X)(cutree(dend, h=X) %>% unique %>% length))
df_hclust <- data.frame(Height=height, N.Clusters=n_clusters)
ggplotly(ggplot(df_hclust, aes(x=Height, y=N.Clusters)) + theme_minimal() +
           geom_line(stat='identity', color='darkgreen'))

O gráfico plotado é reativo, podemos encostar o mouse sobre ele e verificar, manualmente, os valores de \(X\) (altura) e \(Y\) (número de clusters) para cada uma das barras. Podemos observar a existência de um platô quando a altura do algoritmo é igual a \(0.2\). Nesse caso, deixamos de trabalhar com \(200\) variáveis e passamos a operar com \(40\) grupos diferentes.

Iremos então acessar o nosso dicionário de metadados e atribuir a cada uma das variáveis o respectivo identificador do cluster.

opt_cut <- cutree(dend, h=0.2)
df_in_clusters <- data.frame('Indicator.Code'=names(opt_cut), 
                             'Cluster.Index'=as.vector(opt_cut))

df_in_clusters <- df_in_clusters %>% 
  inner_join(df_in_meta, by=(c('Indicator.Code'='Indicator Code')))

## Warning: Column `Indicator.Code`/`Indicator Code` joining factor and
## character vector, coercing into character vector

df_in_clusters %>% head

write.csv(df_in_clusters, './Agrupamento_Variaveis.csv')

Finalmente, podemos observar o arquivo de saída e: * Descrever cada um dos grupos * Interpretar e caracterizar cada cluster * Identificar potenciais correlações espúrias

Isso será feito na próxima seção.

3.4. Interpretação dos Clusters e Identificação de Correlações Espúrias

Uma tabela caracterizando cada um dos clusteres foi criada. Ela será impressa nesta seção para consulta e comparação com as análises que serão feitas daqui em diante. As descrições foram feitas relacionando-se cada grupo de variáveis dentro de um mesmo cluster e interpretando o que cada agrupamento representa.

df_grupos <- read_excel('./GRUPOS.xlsx')
df_grupos

3.5. Avaliação da Qualidade dos Clusters

Podemos avaliar a qualidade dos clusters verificando a distribuição das correlações para agrupamento. Queremos grupos com distribuições de correlações o mais próximas de \(100\%\) que for possível.

get_corr_distrib <- function(cluster_index) {
  
  cluster_table <- df_in_clusters %>% filter(Cluster.Index == cluster_index)
  vars_list <- cluster_table[['Indicator.Code']]
  n_cluster_vars <- length(vars_list)
  
  if (n_cluster_vars == 1) {
    return(NA)
  }
  
  corr_vec <- c()
  for (i in 2:n_cluster_vars) {
    for (j in 1:(i - 1)) {
      corr_vec <- c(corr_vec, cor_mat[vars_list[[i]], vars_list[[j]]] %>% abs)
    }
  }
  return(corr_vec)
}

corr_vec <- c()
cluster_vec <- c()

for (k in (df_in_clusters[['Cluster.Index']] %>% unique)) {
  new_rows <- get_corr_distrib(k)
  corr_vec <- c(corr_vec, new_rows)
  cluster_vec <- c(cluster_vec, rep(k %>% as.character, times=new_rows %>% length))
}

df_corr_per_cluster <- data.frame(Correlation=corr_vec, Cluster.Index=cluster_vec) %>%
  mutate(Cluster.Index=as.numeric(Cluster.Index)) %>%
  inner_join(df_grupos, by=c('Cluster.Index'='Cluster.Index'))

ggplot(df_corr_per_cluster) + 
  geom_density(aes(x=Correlation, fill=Cluster.Cod), alpha=.1) +
  scale_x_continuous(limits=c(.7, 1)) + theme_minimal() +
  labs(x='Correlação', y='Densidade por Cluster',
       title='Distribuição da correlação',
       subtitle='Por pares de variáveis dentro de um mesmo grupo / por grupo')

## Warning: Removed 17 rows containing non-finite values (stat_density).

## Warning: Groups with fewer than two data points have been dropped.

## Warning: Groups with fewer than two data points have been dropped.

## Warning: Groups with fewer than two data points have been dropped.

## Warning: Groups with fewer than two data points have been dropped.

## Warning: Groups with fewer than two data points have been dropped.

Podemos ainda plotar a distribuição global, sem segregação por cluster:

ggplot(df_corr_per_cluster, aes(x=Correlation)) + geom_density() +
  scale_x_continuous(limits=c(.7, 1)) + theme_minimal() +
  labs(x='Correlação (Variáveis em um mesmo cluster)', y='Densidade',
       title='Distribuição da correlação',
       subtitle='Por pares de variáveis dentro de um mesmo grupo / geral')

## Warning: Removed 17 rows containing non-finite values (stat_density).

Podemos então observar que, além de reduzir o número de variáveis de \(200\) para \(32\), garantimos que dentro de um cluster teremos correlações iguais a, no mínimo, \(70\%\). Assim, podemos prosseguir nossa análise, focando apenas em variáveis agregadas por grupos. A agregação será realizada em duas etapas:

• Variáveis com correlação maior ou igual a \(99\%\) podem ser consideradas como sendo a mesma variável. Assim, iremos em uma primeira etapa retirar a média de grupos de variáveis com tal correlação fortíssima.

• Em seguida, agruparemos o cluster realizando, também, uma média aritimética. Nesta etapa, os erros gerados por cada variável também serão submetidos a tal média e isso contribui na normalização do erro (pela lei do limite central), na redução da variância do erro (pois o desvio padrão se reduz em um fator igual a \(\sqrt{N_{Cluster}}\)) e na eliminação do problema da multicolinearidade sem eliminarmos qualquer variável da análise.

• O agrupamento por média deverá ser efetuado sobre as variáveis normalizadas, para que o efeito da média não seja comparado a uma média ponderada na qual variáveis com unidades de maiores ordens de grandeza prevalecem.

df_in_p_sd <- df_in_p
df_in_p_sd <- df_in_p_sd %>%
  group_by(Indicator) %>%
  summarise_at(vars(Value), sd) %>%
  rename(Sd = Value) %>%
  as.data.frame

df_in_p_mean <- df_in_p
df_in_p_mean <- df_in_p_mean %>%
  group_by(Indicator) %>%
  summarise_at(vars(Value), mean) %>%
  rename(Mean = Value) %>%
  as.data.frame

df_in_p_grupos <- df_in_p %>%
  mutate(Indicator = as.character(Indicator)) %>%
  inner_join(df_in_p_sd, by=c('Indicator'='Indicator')) %>%
  inner_join(df_in_p_mean, by=c('Indicator' = 'Indicator')) %>%
  mutate(Value = (Value - Mean) / Sd) %>%
  inner_join(df_in_clusters, by=c('Indicator' = 'Indicator.Code')) %>%
  mutate(Cluster.Index = as.numeric(Cluster.Index)) %>%
  inner_join(df_grupos, by=c('Cluster.Index' = 'Cluster.Index')) %>%
  select('Cluster.Index', 'Cluster.Cod', 'Cluster.Description', 'Year', 'Value') %>%
  group_by(Cluster.Cod, Year) %>%
  summarise_at(vars(Value), mean) %>%
  as.data.frame

df_in_p_grupos %>% head

Há correlação entre grupos de um mesmo cluster? É necessário fazer essa verificação final para conferir se o problema da multicolinearidade pode ocorrer ainda. Iniciemos plotando a evolução do valor médio encontrado em cada um dos clusters:

clusters_list <- df_in_p_grupos[['Cluster.Cod']] %>% unique

Gráficos das Variáveis

Variáveis 1 a 6

ggplot(df_in_p_grupos %>% filter(Cluster.Cod %in% clusters_list[1:6]), aes(x=Year, y=Value, color=Cluster.Cod)) + geom_line() + theme_minimal() +
  labs(x='Ano', y ='Valor', title='Evolução temporal das variáveis',
       subtitle='Variáveis 1 a 6')

Variáveis 7 a 12

ggplot(df_in_p_grupos %>% filter(Cluster.Cod %in% clusters_list[7:12]), aes(x=Year, y=Value, color=Cluster.Cod)) + geom_line() + theme_minimal() +
  labs(x='Ano', y ='Valor', title='Evolução temporal das variáveis',
       subtitle='Variáveis 7 a 12')

Variáveis 13 a 18

ggplot(df_in_p_grupos %>% filter(Cluster.Cod %in% clusters_list[13:18]), aes(x=Year, y=Value, color=Cluster.Cod)) + geom_line() + theme_minimal() +
  labs(x='Ano', y ='Valor', title='Evolução temporal das variáveis',
       subtitle='Variáveis 13 a 18')

Variáveis 19 a 24

ggplot(df_in_p_grupos %>% filter(Cluster.Cod %in% clusters_list[19:24]), aes(x=Year, y=Value, color=Cluster.Cod)) + geom_line() + theme_minimal() +
  labs(x='Ano', y ='Valor', title='Evolução temporal das variáveis',
       subtitle='Variáveis 19 a 24')

Variáveis 25 a 30

ggplot(df_in_p_grupos %>% filter(Cluster.Cod %in% clusters_list[25:30]), aes(x=Year, y=Value, color=Cluster.Cod)) + geom_line() + theme_minimal() +
  labs(x='Ano', y ='Valor', title='Evolução temporal das variáveis',
       subtitle='Variáveis 25 a 30')

Variáveis 31 a 36

ggplot(df_in_p_grupos %>% filter(Cluster.Cod %in% clusters_list[31:36]), aes(x=Year, y=Value, color=Cluster.Cod)) + geom_line() + theme_minimal() +
  labs(x='Ano', y ='Valor', title='Evolução temporal das variáveis',
       subtitle='Variáveis 31 a 36')

Variáveis 37 a 39

ggplot(df_in_p_grupos %>% filter(Cluster.Cod %in% clusters_list[37:39]), aes(x=Year, y=Value, color=Cluster.Cod)) + geom_line() + theme_minimal() +
  labs(x='Ano', y ='Valor', title='Evolução temporal das variáveis',
       subtitle='Variáveis 37 a 39')

Verificação de Correlograma

Para nos certificarmos de maneira exata, verifiquemos o correlograma. Para isso, precisaremos retirar a tabela utilizada nos gráficos anteriores da forma de pivô.

list_years <- df_in_p_grupos[['Year']] %>% unique
list_clusters <- df_in_p_grupos[['Cluster.Cod']] %>% unique

nrows_mat <- length(list_years)
ncols_mat <- length(list_clusters)

df_in_t_grupos <- zeros(nrows_mat, ncols_mat)

rownames(df_in_t_grupos) <- list_years
colnames(df_in_t_grupos) <- list_clusters

for (curr_year in list_years) {
  for (curr_cluster in list_clusters) {
    df_in_t_grupos[[curr_year %>% as.character, curr_cluster]] <-
      (df_in_p_grupos %>% 
       filter(Year == curr_year &
              Cluster.Cod == curr_cluster))$Value[[1]]
  }
}

df_in_t_grupos <- df_in_t_grupos %>% as.data.frame
mat_cor_grupos <- cor(df_in_t_grupos) %>% abs

mask <- zeros(ncols_mat) > 0
for (i in 1:ncols_mat) {
  mat_cor_grupos[[i, i]] <- 0
  mask[[i]] <- (any(abs(mat_cor_grupos[i,]) > 0.78))
  mat_cor_grupos[[i, i]] <- 1
}

mat_cor_grupos[mask, mask] %>% corrplot

Ainda temos um grande problema de multicolinearidade entre os clusters. Para resolvermos isso, podemos recorrer à técnica da análise de componentes principais.

3.6. Solução da Multicolinearidade entre Clusters: Análise de Componentes Principais (PCA)

Poderíamos, desde o início, ter realizado essa técnica, sem o intermédio de qualquer algoritmo de clusterização. No entanto, analisar a influência de um número reduzido de clusters em cada um dos componentes principais é mais fácil e garante uma maior interpretabilidade ao modelo.

pr_comp <- prcomp(df_in_t_grupos)
list_var <- (pr_comp$sdev ** 2) / (sum(pr_comp$sdev ** 2))
list_cum <- cumsum(list_var)
df_pr <- data.frame(Perc.Var=list_var,
                    Cum.Perc=list_cum,
                    PC.Index=1:length(list_var))

ggplot(df_pr, aes(x=PC.Index)) + 
  geom_bar(aes(y=Perc.Var), stat='identity', color='darkgreen', fill='white') + 
  geom_line(aes(y=Cum.Perc)) + 
  scale_x_continuous(breaks = seq(0, 40, len=10)) +
  theme_minimal() +
  labs(x='Índice da Componente Principal',
       y='Percentual de Importância da Componente',
       title='Gráfico de Pareto: Componentes Principais',
       subtitle='Importância relativa (barras) e importância cumulativa (linha)')

pr_comp %>% summary

## Importance of components:
##                           PC1    PC2   PC3     PC4     PC5     PC6    PC7
## Standard deviation     3.1518 2.1739 2.005 1.46340 1.22961 1.10821 1.0430
## Proportion of Variance 0.3361 0.1599 0.136 0.07245 0.05115 0.04155 0.0368
## Cumulative Proportion  0.3361 0.4960 0.632 0.70441 0.75556 0.79711 0.8339
##                           PC8     PC9    PC10    PC11    PC12    PC13
## Standard deviation     0.9650 0.87051 0.79517 0.68265 0.58170 0.56293
## Proportion of Variance 0.0315 0.02564 0.02139 0.01577 0.01145 0.01072
## Cumulative Proportion  0.8654 0.89105 0.91244 0.92821 0.93965 0.95037
##                           PC14   PC15   PC16    PC17    PC18    PC19
## Standard deviation     0.49146 0.4582 0.4384 0.40714 0.38020 0.32840
## Proportion of Variance 0.00817 0.0071 0.0065 0.00561 0.00489 0.00365
## Cumulative Proportion  0.95855 0.9657 0.9721 0.97776 0.98265 0.98630
##                           PC20    PC21    PC22    PC23    PC24    PC25
## Standard deviation     0.28852 0.26693 0.24880 0.20955 0.19867 0.16624
## Proportion of Variance 0.00282 0.00241 0.00209 0.00149 0.00134 0.00093
## Cumulative Proportion  0.98911 0.99152 0.99362 0.99510 0.99644 0.99737
##                           PC26    PC27    PC28    PC29    PC30    PC31
## Standard deviation     0.14316 0.11484 0.10626 0.09295 0.08766 0.06892
## Proportion of Variance 0.00069 0.00045 0.00038 0.00029 0.00026 0.00016
## Cumulative Proportion  0.99807 0.99851 0.99890 0.99919 0.99945 0.99961
##                           PC32    PC33    PC34    PC35    PC36    PC37
## Standard deviation     0.06533 0.05669 0.04332 0.03426 0.02667 0.01283
## Proportion of Variance 0.00014 0.00011 0.00006 0.00004 0.00002 0.00001
## Cumulative Proportion  0.99975 0.99986 0.99992 0.99996 0.99999 0.99999
##                           PC38     PC39
## Standard deviation     0.01123 0.007274
## Proportion of Variance 0.00000 0.000000
## Cumulative Proportion  1.00000 1.000000

Podemos observar que grande parte das variações nas séries pode ser explicada por um pequeno número de componentes principais (conseguimos explicar 95% das variações a partir de 12 componentes). Assim, esperamos que, na etapa de modelagem, consigamos criar um modelo com poucas variáveis explicativas e buscaremos eliminá-las a partir de duas etapas:

• Por meio da função STEP, que buscará eliminar coeficientes de forma a reduzir o parâmetro AIC (Akaike Information Criterion)
• Em seguida, verificaremos quais coeficientes na saída são pouco significantes e possuem altos valores-P nos testes de Student.

Em cada instante de tempo, o vetor \(\mathcal{C} = (C_1, C_2,..., C_{39})\) formado pelos clusteres será multiplicado pela matriz de rotação encontrada na análise de componentes principais e, assim, encontraremos novas séries temporais sobre as quais serão efetuadas as análises de regressão.

Aplicando a transformação em nossa matriz de clusters, obtemos o termo “x” da saída da função prcomp.

df_pca <- pr_comp$x %>% as.data.frame
df_pca$Y <- df_out

Finalmente, podemos iniciar a modelagem do regressor linear.

4. Aplicação e Análise do Modelo de Regressão

4.1. Encontrando Modelo Adequado

Ajustando o modelo, nas componentes principais:

mod = lm(Y~., data=df_pca)
summary(mod)

## 
## Call:
## lm(formula = Y ~ ., data = df_pca)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6328.6 -1902.9  -338.7  1546.3  6104.5 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  391149.2      674.6 579.861  < 2e-16 ***
## PC1            7948.6      216.0  36.807  < 2e-16 ***
## PC2            1689.9      313.1   5.397 5.93e-05 ***
## PC3          -20651.3      339.5 -60.832  < 2e-16 ***
## PC4          -19536.4      465.1 -42.003  < 2e-16 ***
## PC5           -5201.2      553.6  -9.396 6.49e-08 ***
## PC6            8640.3      614.2  14.068 1.99e-10 ***
## PC7             419.2      652.6   0.642 0.529784    
## PC8             733.6      705.3   1.040 0.313762    
## PC9           -2792.0      781.9  -3.571 0.002552 ** 
## PC10           -889.0      856.0  -1.039 0.314453    
## PC11          -7610.6      997.1  -7.633 1.01e-06 ***
## PC12          -4530.4     1170.1  -3.872 0.001352 ** 
## PC13           8178.5     1209.2   6.764 4.55e-06 ***
## PC14          -6912.5     1385.0  -4.991 0.000133 ***
## PC15           6402.1     1485.6   4.310 0.000540 ***
## PC16           2178.5     1552.6   1.403 0.179703    
## PC17          -3261.5     1671.8  -1.951 0.068809 .  
## PC18          -3827.9     1790.3  -2.138 0.048277 *  
## PC19           1670.9     2072.6   0.806 0.431963    
## PC20          -3589.7     2359.1  -1.522 0.147613    
## PC21          15742.9     2550.0   6.174 1.34e-05 ***
## PC22           9454.3     2735.8   3.456 0.003254 ** 
## PC23          15197.6     3248.2   4.679 0.000252 ***
## PC24          16195.3     3426.1   4.727 0.000228 ***
## PC25           6123.6     4094.4   1.496 0.154217    
## PC26          13572.9     4754.5   2.855 0.011468 *  
## PC27           1439.1     5927.0   0.243 0.811247    
## PC28          -6166.5     6405.3  -0.963 0.350020    
## PC29          41051.8     7323.0   5.606 3.94e-05 ***
## PC30          12587.3     7764.7   1.621 0.124535    
## PC31          37646.4     9876.8   3.812 0.001535 ** 
## PC32          -9955.8    10418.3  -0.956 0.353486    
## PC33         -32279.3    12006.6  -2.688 0.016149 *  
## PC34          15625.8    15713.2   0.994 0.334811    
## PC35          20042.1    19866.4   1.009 0.328062    
## PC36          14107.2    25520.1   0.553 0.588049    
## PC37         102290.8    53069.9   1.927 0.071866 .  
## PC38        -108663.4    60589.4  -1.793 0.091824 .  
## PC39         -78836.9    93575.7  -0.842 0.411929    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5048 on 16 degrees of freedom
## Multiple R-squared:  0.9979, Adjusted R-squared:  0.9927 
## F-statistic: 192.3 on 39 and 16 DF,  p-value: 7.613e-16

Observamos que para diversas variáveis não podemos rejeitar a hipótese \(\mathcal{H}_0\) de nulidade do coeficiente. Assim, iremos tentar resolver o problema eliminando as variáveis em sequência adequada a aumentar a verossimilhança do modelo por meio da função “step”, utilizando-se o parâmetro “direction = backward”:

mod_ajustado <- step(mod, direction = 'backward')

## Start:  AIC=964.84
## Y ~ PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC7 + PC8 + PC9 + PC10 + 
##     PC11 + PC12 + PC13 + PC14 + PC15 + PC16 + PC17 + PC18 + PC19 + 
##     PC20 + PC21 + PC22 + PC23 + PC24 + PC25 + PC26 + PC27 + PC28 + 
##     PC29 + PC30 + PC31 + PC32 + PC33 + PC34 + PC35 + PC36 + PC37 + 
##     PC38 + PC39
## 
##        Df  Sum of Sq        RSS     AIC
## - PC27  1 1.5022e+06 4.0921e+08  963.05
## - PC36  1 7.7865e+06 4.1549e+08  963.90
## - PC7   1 1.0512e+07 4.1822e+08  964.26
## <none>               4.0770e+08  964.84
## - PC19  1 1.6561e+07 4.2427e+08  965.07
## - PC39  1 1.8087e+07 4.2579e+08  965.27
## - PC32  1 2.3269e+07 4.3097e+08  965.95
## - PC28  1 2.3616e+07 4.3132e+08  965.99
## - PC34  1 2.5199e+07 4.3290e+08  966.20
## - PC35  1 2.5934e+07 4.3364e+08  966.29
## - PC10  1 2.7484e+07 4.3519e+08  966.49
## - PC8   1 2.7565e+07 4.3527e+08  966.50
## - PC16  1 5.0163e+07 4.5787e+08  969.34
## - PC25  1 5.6999e+07 4.6470e+08  970.17
## - PC20  1 5.9001e+07 4.6671e+08  970.41
## - PC30  1 6.6964e+07 4.7467e+08  971.36
## - PC38  1 8.1959e+07 4.8966e+08  973.10
## - PC37  1 9.4668e+07 5.0237e+08  974.53
## - PC17  1 9.6981e+07 5.0469e+08  974.79
## - PC18  1 1.1650e+08 5.2420e+08  976.91
## - PC33  1 1.8418e+08 5.9188e+08  983.71
## - PC26  1 2.0766e+08 6.1536e+08  985.89
## - PC22  1 3.0431e+08 7.1201e+08  994.06
## - PC9   1 3.2490e+08 7.3261e+08  995.66
## - PC31  1 3.7020e+08 7.7790e+08  999.02
## - PC12  1 3.8197e+08 7.8968e+08  999.86
## - PC15  1 4.7324e+08 8.8094e+08 1005.98
## - PC23  1 5.5779e+08 9.6550e+08 1011.12
## - PC24  1 5.6938e+08 9.7708e+08 1011.78
## - PC14  1 6.3476e+08 1.0425e+09 1015.41
## - PC2   1 7.4226e+08 1.1500e+09 1020.91
## - PC29  1 8.0078e+08 1.2085e+09 1023.69
## - PC21  1 9.7120e+08 1.3789e+09 1031.08
## - PC13  1 1.1658e+09 1.5735e+09 1038.47
## - PC11  1 1.4845e+09 1.8922e+09 1048.80
## - PC5   1 2.2496e+09 2.6573e+09 1067.81
## - PC6   1 5.0428e+09 5.4505e+09 1108.04
## - PC1   1 3.4521e+10 3.4928e+10 1212.07
## - PC4   1 4.4955e+10 4.5363e+10 1226.71
## - PC3   1 9.4295e+10 9.4702e+10 1267.92
## 
## Step:  AIC=963.05
## Y ~ PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC7 + PC8 + PC9 + PC10 + 
##     PC11 + PC12 + PC13 + PC14 + PC15 + PC16 + PC17 + PC18 + PC19 + 
##     PC20 + PC21 + PC22 + PC23 + PC24 + PC25 + PC26 + PC28 + PC29 + 
##     PC30 + PC31 + PC32 + PC33 + PC34 + PC35 + PC36 + PC37 + PC38 + 
##     PC39
## 
##        Df  Sum of Sq        RSS     AIC
## - PC36  1 7.7865e+06 4.1699e+08  962.10
## - PC7   1 1.0512e+07 4.1972e+08  962.47
## <none>               4.0921e+08  963.05
## - PC19  1 1.6561e+07 4.2577e+08  963.27
## - PC39  1 1.8087e+07 4.2729e+08  963.47
## - PC32  1 2.3269e+07 4.3248e+08  964.14
## - PC28  1 2.3616e+07 4.3282e+08  964.19
## - PC34  1 2.5199e+07 4.3441e+08  964.39
## - PC35  1 2.5934e+07 4.3514e+08  964.49
## - PC10  1 2.7484e+07 4.3669e+08  964.69
## - PC8   1 2.7565e+07 4.3677e+08  964.70
## - PC16  1 5.0163e+07 4.5937e+08  967.52
## - PC25  1 5.6999e+07 4.6621e+08  968.35
## - PC20  1 5.9001e+07 4.6821e+08  968.59
## - PC30  1 6.6964e+07 4.7617e+08  969.53
## - PC38  1 8.1959e+07 4.9117e+08  971.27
## - PC37  1 9.4668e+07 5.0387e+08  972.70
## - PC17  1 9.6981e+07 5.0619e+08  972.96
## - PC18  1 1.1650e+08 5.2570e+08  975.07
## - PC33  1 1.8418e+08 5.9338e+08  981.86
## - PC26  1 2.0766e+08 6.1687e+08  984.03
## - PC22  1 3.0431e+08 7.1351e+08  992.18
## - PC9   1 3.2490e+08 7.3411e+08  993.77
## - PC31  1 3.7020e+08 7.7941e+08  997.13
## - PC12  1 3.8197e+08 7.9118e+08  997.97
## - PC15  1 4.7324e+08 8.8245e+08 1004.08
## - PC23  1 5.5779e+08 9.6700e+08 1009.20
## - PC24  1 5.6938e+08 9.7858e+08 1009.87
## - PC14  1 6.3476e+08 1.0440e+09 1013.49
## - PC2   1 7.4226e+08 1.1515e+09 1018.98
## - PC29  1 8.0078e+08 1.2100e+09 1021.76
## - PC21  1 9.7120e+08 1.3804e+09 1029.14
## - PC13  1 1.1658e+09 1.5750e+09 1036.52
## - PC11  1 1.4845e+09 1.8937e+09 1046.84
## - PC5   1 2.2496e+09 2.6588e+09 1065.84
## - PC6   1 5.0428e+09 5.4520e+09 1106.06
## - PC1   1 3.4521e+10 3.4930e+10 1210.07
## - PC4   1 4.4955e+10 4.5364e+10 1224.71
## - PC3   1 9.4295e+10 9.4704e+10 1265.93
## 
## Step:  AIC=962.1
## Y ~ PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC7 + PC8 + PC9 + PC10 + 
##     PC11 + PC12 + PC13 + PC14 + PC15 + PC16 + PC17 + PC18 + PC19 + 
##     PC20 + PC21 + PC22 + PC23 + PC24 + PC25 + PC26 + PC28 + PC29 + 
##     PC30 + PC31 + PC32 + PC33 + PC34 + PC35 + PC37 + PC38 + PC39
## 
##        Df  Sum of Sq        RSS     AIC
## - PC7   1 1.0512e+07 4.2751e+08  961.50
## <none>               4.1699e+08  962.10
## - PC19  1 1.6561e+07 4.3355e+08  962.28
## - PC39  1 1.8087e+07 4.3508e+08  962.48
## - PC32  1 2.3269e+07 4.4026e+08  963.14
## - PC28  1 2.3616e+07 4.4061e+08  963.19
## - PC34  1 2.5199e+07 4.4219e+08  963.39
## - PC35  1 2.5934e+07 4.4293e+08  963.48
## - PC10  1 2.7484e+07 4.4448e+08  963.68
## - PC8   1 2.7565e+07 4.4456e+08  963.69
## - PC16  1 5.0163e+07 4.6716e+08  966.46
## - PC25  1 5.6999e+07 4.7399e+08  967.28
## - PC20  1 5.9001e+07 4.7599e+08  967.51
## - PC30  1 6.6964e+07 4.8396e+08  968.44
## - PC38  1 8.1959e+07 4.9895e+08  970.15
## - PC37  1 9.4668e+07 5.1166e+08  971.56
## - PC17  1 9.6981e+07 5.1397e+08  971.81
## - PC18  1 1.1650e+08 5.3349e+08  973.90
## - PC33  1 1.8418e+08 6.0117e+08  980.59
## - PC26  1 2.0766e+08 6.2465e+08  982.73
## - PC22  1 3.0431e+08 7.2130e+08  990.79
## - PC9   1 3.2490e+08 7.4190e+08  992.36
## - PC31  1 3.7020e+08 7.8719e+08  995.68
## - PC12  1 3.8197e+08 7.9896e+08  996.51
## - PC15  1 4.7324e+08 8.9023e+08 1002.57
## - PC23  1 5.5779e+08 9.7479e+08 1007.65
## - PC24  1 5.6938e+08 9.8637e+08 1008.31
## - PC14  1 6.3476e+08 1.0517e+09 1011.91
## - PC2   1 7.4226e+08 1.1593e+09 1017.36
## - PC29  1 8.0078e+08 1.2178e+09 1020.12
## - PC21  1 9.7120e+08 1.3882e+09 1027.45
## - PC13  1 1.1658e+09 1.5828e+09 1034.80
## - PC11  1 1.4845e+09 1.9015e+09 1045.07
## - PC5   1 2.2496e+09 2.6666e+09 1064.01
## - PC6   1 5.0428e+09 5.4597e+09 1104.14
## - PC1   1 3.4521e+10 3.4938e+10 1208.08
## - PC4   1 4.4955e+10 4.5372e+10 1222.72
## - PC3   1 9.4295e+10 9.4712e+10 1263.93
## 
## Step:  AIC=961.5
## Y ~ PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC8 + PC9 + PC10 + PC11 + 
##     PC12 + PC13 + PC14 + PC15 + PC16 + PC17 + PC18 + PC19 + PC20 + 
##     PC21 + PC22 + PC23 + PC24 + PC25 + PC26 + PC28 + PC29 + PC30 + 
##     PC31 + PC32 + PC33 + PC34 + PC35 + PC37 + PC38 + PC39
## 
##        Df  Sum of Sq        RSS     AIC
## <none>               4.2751e+08  961.50
## - PC19  1 1.6561e+07 4.4407e+08  961.62
## - PC39  1 1.8087e+07 4.4559e+08  961.82
## - PC32  1 2.3269e+07 4.5077e+08  962.46
## - PC28  1 2.3616e+07 4.5112e+08  962.51
## - PC34  1 2.5199e+07 4.5270e+08  962.70
## - PC35  1 2.5934e+07 4.5344e+08  962.79
## - PC10  1 2.7484e+07 4.5499e+08  962.98
## - PC8   1 2.7565e+07 4.5507e+08  962.99
## - PC16  1 5.0163e+07 4.7767e+08  965.71
## - PC25  1 5.6999e+07 4.8450e+08  966.50
## - PC20  1 5.9001e+07 4.8651e+08  966.73
## - PC30  1 6.6964e+07 4.9447e+08  967.64
## - PC38  1 8.1959e+07 5.0946e+08  969.32
## - PC37  1 9.4668e+07 5.2217e+08  970.70
## - PC17  1 9.6981e+07 5.2449e+08  970.94
## - PC18  1 1.1650e+08 5.4400e+08  972.99
## - PC33  1 1.8418e+08 6.1168e+08  979.56
## - PC26  1 2.0766e+08 6.3516e+08  981.67
## - PC22  1 3.0431e+08 7.3181e+08  989.60
## - PC9   1 3.2490e+08 7.5241e+08  991.15
## - PC31  1 3.7020e+08 7.9771e+08  994.43
## - PC12  1 3.8197e+08 8.0948e+08  995.25
## - PC15  1 4.7324e+08 9.0074e+08 1001.23
## - PC23  1 5.5779e+08 9.8530e+08 1006.25
## - PC24  1 5.6938e+08 9.9688e+08 1006.91
## - PC14  1 6.3476e+08 1.0623e+09 1010.47
## - PC2   1 7.4226e+08 1.1698e+09 1015.86
## - PC29  1 8.0078e+08 1.2283e+09 1018.60
## - PC21  1 9.7120e+08 1.3987e+09 1025.87
## - PC13  1 1.1658e+09 1.5933e+09 1033.17
## - PC11  1 1.4845e+09 1.9120e+09 1043.38
## - PC5   1 2.2496e+09 2.6771e+09 1062.23
## - PC6   1 5.0428e+09 5.4703e+09 1102.25
## - PC1   1 3.4521e+10 3.4948e+10 1206.10
## - PC4   1 4.4955e+10 4.5383e+10 1220.73
## - PC3   1 9.4295e+10 9.4722e+10 1261.94

summary(mod_ajustado)

## 
## Call:
## lm(formula = Y ~ PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC8 + PC9 + 
##     PC10 + PC11 + PC12 + PC13 + PC14 + PC15 + PC16 + PC17 + PC18 + 
##     PC19 + PC20 + PC21 + PC22 + PC23 + PC24 + PC25 + PC26 + PC28 + 
##     PC29 + PC30 + PC31 + PC32 + PC33 + PC34 + PC35 + PC37 + PC38 + 
##     PC39, data = df_pca)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6119.0 -1972.8  -442.5  1677.2  6079.0 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  391149.2      633.9 617.082  < 2e-16 ***
## PC1            7948.6      202.9  39.169  < 2e-16 ***
## PC2            1689.9      294.2   5.744 1.55e-05 ***
## PC3          -20651.3      319.0 -64.737  < 2e-16 ***
## PC4          -19536.4      437.1 -44.699  < 2e-16 ***
## PC5           -5201.2      520.2  -9.999 5.27e-09 ***
## PC6            8640.3      577.2  14.971 5.70e-12 ***
## PC8             733.6      662.8   1.107 0.282176    
## PC9           -2792.0      734.7  -3.800 0.001210 ** 
## PC10           -889.0      804.4  -1.105 0.282864    
## PC11          -7610.6      936.9  -8.123 1.34e-07 ***
## PC12          -4530.4     1099.5  -4.120 0.000582 ***
## PC13           8178.5     1136.2   7.198 7.75e-07 ***
## PC14          -6912.5     1301.4  -5.311 3.98e-05 ***
## PC15           6402.1     1396.0   4.586 0.000202 ***
## PC16           2178.5     1459.0   1.493 0.151822    
## PC17          -3261.5     1571.0  -2.076 0.051699 .  
## PC18          -3827.9     1682.3  -2.275 0.034651 *  
## PC19           1670.9     1947.6   0.858 0.401630    
## PC20          -3589.7     2216.8  -1.619 0.121857    
## PC21          15742.9     2396.2   6.570 2.73e-06 ***
## PC22           9454.3     2570.8   3.678 0.001599 ** 
## PC23          15197.6     3052.3   4.979 8.33e-05 ***
## PC24          16195.3     3219.5   5.030 7.43e-05 ***
## PC25           6123.6     3847.4   1.592 0.127971    
## PC26          13572.9     4467.8   3.038 0.006768 ** 
## PC28          -6166.5     6019.0  -1.025 0.318466    
## PC29          41051.8     6881.3   5.966 9.66e-06 ***
## PC30          12587.3     7296.4   1.725 0.100730    
## PC31          37646.4     9281.1   4.056 0.000674 ***
## PC32          -9955.8     9789.9  -1.017 0.321953    
## PC33         -32279.3    11282.4  -2.861 0.009998 ** 
## PC34          15625.8    14765.4   1.058 0.303201    
## PC35          20042.1    18668.1   1.074 0.296447    
## PC37         102290.8    49868.9   2.051 0.054295 .  
## PC38        -108663.4    56934.8  -1.909 0.071543 .  
## PC39         -78836.9    87931.4  -0.897 0.381166    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4743 on 19 degrees of freedom
## Multiple R-squared:  0.9978, Adjusted R-squared:  0.9935 
## F-statistic: 235.9 on 36 and 19 DF,  p-value: < 2.2e-16

Encontramos um modelo com apenas \(36\) variáveis, grande parte delas é significante e apenas uma eliminação foi realizada. O coeficiente \(R^2\) é altíssimo, na ordem de \(99%\), indicando, a princípio, que uma grande parcela da variação quadrática total pode, de fato, ser explicada pelo modelo.

Porém, para termos uma noção geral dos resultados, precisamos olhar mais detalhadamente os resíduos.

4.2. Pré-Análise de Resíduos

Observemos os resíduos do modelo ajustado:

par(mfrow=c(2, 2))
mod_ajustado %>% aov %>% plot

• A curva de resíduos e desvios-padrões por valores ajustados não demonstra nenhum padrão, o que indica normalidade e ausência de autocorrelação.
• A curva quantil-quantil gaussiana forma um padrão de reta de \(45\) graus, o que se encontra de acordo com a literatura em amostras de erros que seguem distribuições normais
• A última curva, demonstra a inexistência de outliers significativos (que possuam elevadas distâncias de Cook). Assim, a retirada individual de cada ponto da amostra não pode interferir significativamente no modelo já encontrado.

Testes efetuados

anares <- resid(mod_ajustado)
ad.test(anares)

## 
##  Anderson-Darling normality test
## 
## data:  anares
## A = 0.38975, p-value = 0.3719

shapiro.test(anares)

## 
##  Shapiro-Wilk normality test
## 
## data:  anares
## W = 0.98017, p-value = 0.4828

Os p-valores são muito maiores que \(0.05\) (tanto para o teste de Anderson-Darling quanto para o teste de Shapiro-Wilk) e, assim, podemos adotar a hipótese de normalidade. Testando a homocedasticidade por meio do teste de Breusch-Pagan

bptest(mod_ajustado)

## 
##  studentized Breusch-Pagan test
## 
## data:  mod_ajustado
## BP = 34.627, df = 36, p-value = 0.5339

O elevado P-Valor não nos permite rejeitar a hipótese de homocedasticidade. Finalmente, testemos a autocorrelação por meio de um teste de Durblin-Watson:

dwtest(mod_ajustado)

## 
##  Durbin-Watson test
## 
## data:  mod_ajustado
## DW = 2.2876, p-value = 0.1314
## alternative hypothesis: true autocorrelation is greater than 0

Como o p-valor é maior que \(5 \%\). Assim, temos todas as hipóteses necessárias para uma boa regressão linear (já sabemos que não há multicolinearidade pois realizamos decomposição por componentes principais na transformação dos dados).

Nosso modelo cumpre os requisitos de uma boa regressão linear e podemos, agora, interpretar os resultados. Porém, para isso, precisamos verificar a influência de cada um dos clusters sobre cada um dos componentes principais. Isso será realizado na próxima seção.

5. Interpretando o Modelo

Para interpretarmos o modelo, iremos, primeiramente, relembrar a sequência de passos adotada até aqui:

1. Inspeção das variáveis
2. Agrupamento das variáveis em clusteres
3. Criação de novas variáveis: cada cluster é transformado em uma única variável, representada pela média dos membros do agrupamento, cada um deles é normalizado (dividido pelo desvio-padrão) antes da operação.
4. Transformação das variáveis utilizando análise de componentes principais (PCA) para eliminar o problema de multicolinearidade entre os clusters.
5. Aplicação do modelo linear

Assim, cada componente principal será explicado, em diferentes graus e percentuais, por diferentes clusters:

pr_comp_contrib <- pr_comp$rotation ** 2 %>% 
  apply(MARGIN=1, FUN=function(X) (X / sum(X))) %>% t %>% as.data.frame

pr_comp_contrib <- rownames_to_column(pr_comp_contrib, var='Cluster.Name')

Iremos demonstrar as contribuições para os valores principais mais significativos do modelo final:

summary(mod_ajustado)

## 
## Call:
## lm(formula = Y ~ PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC8 + PC9 + 
##     PC10 + PC11 + PC12 + PC13 + PC14 + PC15 + PC16 + PC17 + PC18 + 
##     PC19 + PC20 + PC21 + PC22 + PC23 + PC24 + PC25 + PC26 + PC28 + 
##     PC29 + PC30 + PC31 + PC32 + PC33 + PC34 + PC35 + PC37 + PC38 + 
##     PC39, data = df_pca)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6119.0 -1972.8  -442.5  1677.2  6079.0 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  391149.2      633.9 617.082  < 2e-16 ***
## PC1            7948.6      202.9  39.169  < 2e-16 ***
## PC2            1689.9      294.2   5.744 1.55e-05 ***
## PC3          -20651.3      319.0 -64.737  < 2e-16 ***
## PC4          -19536.4      437.1 -44.699  < 2e-16 ***
## PC5           -5201.2      520.2  -9.999 5.27e-09 ***
## PC6            8640.3      577.2  14.971 5.70e-12 ***
## PC8             733.6      662.8   1.107 0.282176    
## PC9           -2792.0      734.7  -3.800 0.001210 ** 
## PC10           -889.0      804.4  -1.105 0.282864    
## PC11          -7610.6      936.9  -8.123 1.34e-07 ***
## PC12          -4530.4     1099.5  -4.120 0.000582 ***
## PC13           8178.5     1136.2   7.198 7.75e-07 ***
## PC14          -6912.5     1301.4  -5.311 3.98e-05 ***
## PC15           6402.1     1396.0   4.586 0.000202 ***
## PC16           2178.5     1459.0   1.493 0.151822    
## PC17          -3261.5     1571.0  -2.076 0.051699 .  
## PC18          -3827.9     1682.3  -2.275 0.034651 *  
## PC19           1670.9     1947.6   0.858 0.401630    
## PC20          -3589.7     2216.8  -1.619 0.121857    
## PC21          15742.9     2396.2   6.570 2.73e-06 ***
## PC22           9454.3     2570.8   3.678 0.001599 ** 
## PC23          15197.6     3052.3   4.979 8.33e-05 ***
## PC24          16195.3     3219.5   5.030 7.43e-05 ***
## PC25           6123.6     3847.4   1.592 0.127971    
## PC26          13572.9     4467.8   3.038 0.006768 ** 
## PC28          -6166.5     6019.0  -1.025 0.318466    
## PC29          41051.8     6881.3   5.966 9.66e-06 ***
## PC30          12587.3     7296.4   1.725 0.100730    
## PC31          37646.4     9281.1   4.056 0.000674 ***
## PC32          -9955.8     9789.9  -1.017 0.321953    
## PC33         -32279.3    11282.4  -2.861 0.009998 ** 
## PC34          15625.8    14765.4   1.058 0.303201    
## PC35          20042.1    18668.1   1.074 0.296447    
## PC37         102290.8    49868.9   2.051 0.054295 .  
## PC38        -108663.4    56934.8  -1.909 0.071543 .  
## PC39         -78836.9    87931.4  -0.897 0.381166    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4743 on 19 degrees of freedom
## Multiple R-squared:  0.9978, Adjusted R-squared:  0.9935 
## F-statistic: 235.9 on 36 and 19 DF,  p-value: < 2.2e-16

São as componentes: PC1, PC3 e PC4, que possuem significância abaixo de \(0.1\%\) no teste de Student. Plotando as contribuições:

pr_comp_contrib_p <- pr_comp_contrib %>% 
  gather(key='Component', value='Contrib', -Cluster.Name) %>%
  arrange(desc(Contrib))

plot.Radar <- function(component_name) {
  
  options(warn=-1)
  df_to_plot <- pr_comp_contrib_p %>% filter(Component %in% c(component_name))
  
  max_prob <- max(df_to_plot$Contrib)
  size_x <- length(df_to_plot$Cluster.Name %>% unique)
  seq_angle <- 360/(2*pi) * rev(pi/2 + seq(pi/size_x, 
                                           2*pi - pi/14, 
                                           len = size_x))
  
  (ggplot(data=df_to_plot, aes(x=reorder(Cluster.Name, -Contrib), 
                               y=Contrib, group=Component,
                              fill=Contrib)) + 
  geom_point(size=1) + 
  geom_bar(stat='identity') + 
  ggtitle('Análise de Componentes Principais' %>% paste(component_name, sep=': '))  + 
  geom_hline(aes(yintercept=0), lwd=1, lty=2) + 
  scale_y_continuous(limits=c(0, 1.35*max_prob)) + 
  coord_polar() +
  theme_minimal() +
  theme(axis.ticks =element_blank(), 
        axis.text.y =element_blank(), 
        axis.title=element_blank(), 
        axis.text.x=element_text(size = 6, angle=seq_angle))) %>%
    return
}

Composição das Componentes Principais Mais Significativas

Componente 1

plot.Radar('PC1')

Componente 3

plot.Radar('PC3')

Componente 4

plot.Radar('PC4')

Interpretação

Interpretando-se as componentes principais, obtemos:

• A componente principal mais importante é PC1 (maior percentual de variância) possui participações de todas as variáveis mas possui como variáveis mais influentes os clusters: (1) EXIMBEMS (Exportação de Mercadorias, Bens e Serviços), (2) IMPODCEX (Mercadorias Importadas e Alimentos Exportados) e (3) CO2ELAQP (Emissão de CO2 para produção de eletricidade e produção).
• A segunda componente principal analisada, PC3, possui duas variáveis relevantes em sua composição: EXPARMS e CRESPRF, que são, respectivamente, exportação de armas e crescimento populacional em razão da fertilidade.
• A terceira componente principal analisada, PC4, possui como \(3\) variáveis mais relevantes, a Exportação de Bens, Mercadorias e Serviços (EXPBEMS - sem contar com a importação, como no caso da primeira componente), IMPBEMS - a importação de bens, mercadorias e serviços e MIMPOR - Mercadorias Importadas para o Desenvolvimento da América Latina.

Podemos então notar que o comércio internacional possui grande peso na determinação da emissão de CO2 pela França e, em expecial, o comércio de bens importados da América Latina.

Além disso, observamos uma participação relevante da produção de energia e da indústria bélica em tal variável. Logicamente, todos esses fatores estão correlacionados e representam faces diferentes de tal problema ambiental.

6. Interpretando com Wordclouds

A análise das componentes principais na seção anterior nos forneceu alguns insights mas ainda não pareceu ser o suficiente para interpretarmos todas as variáveis explicativas (que são muitas) simultaneamente.

Por isso, propomos um método alternativo para mapear as variáveis relevantes do modelo: iremos gerar um wordcloud com as palavras com maior peso na determinação da saída (emissão de CO2 em KT).

Iremos criar um dataframe relacionando cada palavra contida na descrição das variáveis ao peso de tal variável na determinação da variável explicada. E como determinaremos tal peso?

Iremos utilizar a regra da cadeia para calcular a sensibilidade de cada variável na saída.

• Se \(Y\) é a saída do modelo
• \(X_k\) é alguma variável de entrada do modelo (componente principal), representado pelo coeficiente \(A_k\) então:

\(X_k = \sum_{i=0}^{i=N_{Clusters}} P_i.C_i\) e

\(C_k = \frac{\sum_{i=0}^{i=N_{NV_k}} V_{k,i}}{NV_k}\)

Onde:

• \(N_{Clusters}\) é o número de clusters utilizado
• \(P_i\) é o peso em \(\%\) da variância do clusters \(i\) na variável \(X_k\)
• \(C_k\) é o valor do agrupamento por média do cluster \(k\), que possui \(NV_k\) variáveis e
• \(V_{k,i}\) é o valor da \(i\)-ésima variável do cluster \(k\)

Nesse caso temos:

\(\frac{\partial Y}{\partial V_{k, i}} = \frac{A_k \times P_i}{NV_k}\)

E esse será o peso de cada uma das variáveis, que será somado ao peso de cada uma das palavras contidas dentro da descrição da variável.

Com esses dados, iremos dispor os gráficos em uma nuvem de palavras (iremos utilizar apenas variáveis com p-valor abaixo de \(5 \%\) nessa análise).

concat_descript <- paste(df_in_meta$`Indicator Name`, collapse=' ') %>% toupper()
concat_descript <- gsub('[^[:alnum:] ]', '', concat_descript)

words_vec <- strsplit(concat_descript, split=' ') %>% unlist() %>% unique()
words_weight_vec <- zeros(length(words_vec)) %>% as.vector()
df_wordcloud <- data.frame(Word=words_vec,
                           Weight=words_weight_vec,
                           stringsAsFactors=F)

pc_names <- summary(mod_ajustado)$coefficients %>% rownames()
df_coef_sig <- summary(mod_ajustado)$coefficients %>% as.data.frame()
df_coef_sig$PC <- pc_names
df_coef_sig <- df_coef_sig %>% filter(`Pr(>|t|)` <= 0.05)

df_in_meta_grupos <- df_in_clusters %>% 
  inner_join(df_grupos, by=c('Cluster.Index'='Cluster.Index')) %>%
  inner_join(df_in_meta, by=c('Indicator.Code'='Indicator Code'))

for (var_X in df_coef_sig$PC[-1]) {
  var_A <- mod_ajustado$coefficients[[var_X]]
  for (var_V in pr_comp$rotation[, var_X] %>% names()) {
    
    weight_V <- pr_comp$rotation[[var_V, var_X]] * var_A
    
    V_description <- (df_in_meta_grupos %>% 
                        filter(Cluster.Cod==var_V) %>% 
                        as.data.frame())$`Indicator Name.x`[[1]]
    
    V_description <- gsub('[^[:alnum:] ]', '', V_description)
    
    V_description_words <- V_description %>% toupper() %>% strsplit(split = ' ')
    
    for (word_from_V in V_description_words %>% unlist) {
      curr_val <- (df_wordcloud %>% filter(Word == word_from_V))$Weight[[1]]
      setDT(df_wordcloud)[Word == word_from_V, Weight := curr_val + weight_V]
    }
  }
}

ignore_list <- c('', 'OF', 'PER', 'ANNUAL', 'NET', 
                 'AND', 'OF', 'LCU', 'FROM', 'CONSTANT',
                 'TERMS', 'ADJUSTMENT', 'ON', 'TO', 'KM', 'KWH',
                 'SQ', 'DEC', 'IN', 'TOTAL', 'FINAL', 'BY')

df_wordcloud <- df_wordcloud %>% 
  filter(Weight != 0 & ! Word %in% ignore_list) %>% 
  mutate(Weight = abs(Weight)) %>%
  mutate(Weight=scales::rescale(Weight, to=c(0, 20))) %>%
  arrange(desc(Weight))


df_wordcloud

wordcloud(df_wordcloud$Word, 
          df_wordcloud$Weight,
          colors = brewer.pal(8, 'Dark2'))

• Óleo, produção, gás e combustíveis parecem exercer um primeiro plano nas palavras-chave que possuem mais peso na explicação das emissões de CO2.

• População, bens e trocas também exercem um papel intermediário. Com relação a trocas, temos que, na seção anterior, encontramos relações fortes com diversas variáveis ligadas ao comércio internacional (importações e exportações).

Assim, nosso modelo, ao buscar as variáveis mais relevantes, encontrou três pilares determinantes na explicação das emissões de CO2:

1. Em primeiro lugar: dados relativos à produção energética, principalmente no que diz respeito a Óleo.
2. Em segundo lugar: dados macroeconômicos - população, pessoas, PIB etc.
3. Em terceiro lugar: dados relativos ao comércio internacional (importações e exportações).

E podemos concluir que a emissão de CO2 é um problema desafiante pois se encontra altamente correlacionada com termos-chave relacionados a desenvolvimento e crescimento econômico.

É por isso que a produção precisa ser contida e equilibrada com uma boa diretriz ambiental no conceito de desenvolvimento sustentável.

7. Predição - CO2 em 2014

Para finalizar, iremos prever a emissão de CO2 para o ano de \(2014\) (último ano da análise) conforme o proposto no roteiro deste projeto. Para isso, iremos simplesmente utilizar a função predict.

Iremos também estimar as emissões para os anos de \(2012\) e \(2013\) pois também temos dados perdidos para a emissão de CO2 nesses anos:

df_predict <- predict.lm(object=mod_ajustado,
                       newdata=df_pca[c('2012', '2013', '2014'),],
                       interval='prediction') %>% as.data.frame()

df_predict

Plotando tais resultados:

df_pca$Year <-rownames(df_pca)
df_predict$Year <- c(2012, 2013, 2014)
df_pca_non_missing <- df_pca %>% mutate(Year = as.numeric(Year))
df_pca_non_missing <- df_pca_non_missing %>% filter(Year <= 2012 & Year > 2005)

ggplot() + 
  
  geom_line(data=df_pca_non_missing, 
            aes(x=Year, y=Y, group=1),
            color='blue', stat='identity', group=1) + 
  
  geom_line(data=df_predict, aes(x=Year, y=fit),
            colour='red',
            linetype='dashed') +
  
  geom_linerange(data=df_predict, aes(x=Year, y=fit, ymin=lwr, ymax=upr),
            color='red', group=1) +
  
  geom_point(data=df_pca_non_missing, aes(x=Year, y=Y), 
             shape=21, 
             colour='blue', 
             fill='white', 
             size=3, 
             stroke=1) +
  
  geom_point(data=df_predict, aes(x=Year, y=fit), 
             shape=21, 
             colour='red', 
             fill='white', 
             size=3, 
             stroke=1) +
  
  geom_point(data=df_predict, aes(x=Year, y=lwr), 
             shape=45, 
             colour='red',
             size=15) +
  
  geom_point(data=df_predict, aes(x=Year, y=upr), 
             shape=45, 
             colour='red', 
             size=15) +
  
  labs(
    title = 'Previsão para 2012, 2013 e 2014',
    subtitle = 'Dados históricos em azul, previsões em vermelho',
    x = 'Ano',
    y = 'CO2 (KT) com Erro a 5% (Para previsões)'
  ) + theme_minimal()

O valor esperado para a emissão em KT no ano de \(2014\) é de \(339389.8\) em um intervalo compreendido entre \(327701.9\) e \(351077.6\) e a uma significância de \(5\%\).

Podemos ainda conferir tal resultado com o valor que, de fato ocorreu, para a emissão de CO2 na França nesse ano. Para isso, conferimos esse site: https://www.worldometers.info/co2-emissions/france-co2-emissions/, no qual encontramos uma emissão da ordem de \(320.703.500,00\).

O valor da emissão em \(2014\) se encontra muito acima daquilo que realmente ocorreu e a maior queda na emissão de gás carbônico na França ocorreu justamente nesse ano (\(\approx -9.0\%\)).

É plenamente possível que isso tenha ocorrido devido à implementação de políticas de redução na emissão de gás carbônico que não se encontram diretamente relacionadas às variáveis explicativas do modelo (porém, tal análise foge do escopo do presente projeto).

Em todo caso, isso demonstra que o modelo não é absoluto: os governos podem buscar uma crescente redução em tais emissões no intuito de implementar menos danos ambientais e isso pode ocorrer inclusive em variáveis não incluídas na base de dados aqui utilizada.

Os resultados são positivos para o caso do governo francês, observa-se uma tendêcia de queda.

8. Conclusões

Com base em cada uma das etapas, podemos tirar diversas conclusões:

• No processo de modelagem matemática de determinado problema, a etapa de organização e exploração de dados tendem a tomar uma parcela significativa dos esforços. Essas duas etapas podem ser realizadas conjuntamente pois a organização ocorre ao longo da exploração dos dados brutos que devem ser, paulatinamente, refinados.
• A clusterizaçao tomou um papel importante na etapa de análise exploratória e organização de dados descrito no ponto anterior pois foi fundamental no estudo da redução do número de variáveis perdas significativas de informações.
• Todos os testes de hipóteses foram positivos e os requisitos do modelo foram cumpridos.
• O procedimento de decomposição por valores singulares permitiu, de forma direta, a eliminação do problema de multicolinearidade (que poderia comprometer os resultados de nossos testes de hipóteses). Porém, a aplicação de tal modelo tornou a interpretação do modelo um pouco mais trabalhosa. Entretanto, foi possível verificar a participação das variáveis iniciais nas componentes principais sem maiores prejuízos na interpretabilidade do modelo.
• Emissões de CO2 são altamente relacionadas à indústria de Óleo e Gás e a palavras-chave ligadas à produção energética. Outros aspectos relevantes foram identificados como variáveis ligadas ao comércio mundial, à indústria bélica ou a dados demográficos e macroeconômicos.